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Abstract 

We study the dynamics of a liquid droplet inside a gas over a large range 
of the Knudsen numbers. The moving liquid droplet is modeled by the in- 
compressible Navier-Stokes equations, the surrounding rarefied gas by the 
Boltzmann equation. The interface boundary conditions between the gas 
and liquid phases are derived. The incompressible Navier-Stokes equations 
are solved by a meshfree Lagrangian particle method called Finite Pointset 
Method (FPM), and the Boltzmann equation by a DSMC type of particle 
method. To validiate the coupled solutions of the Boltzmann and the incom- 
pressible Navier-Stokes equations we have further solved the compressible 
and the incompressible Navier-Stokes equations in the gas and liquid phases, 
respectively. In the latter case both the compressible and the incompress- 
ible Navier-Stokes equations are also solved by the FPM. In the continuum 
regime the coupled solutions obtained from the Boltzmann and the incom- 
pressible Navier-Stokes equations match with the solutions obtained from the 
compressible and the incompressible Navier-Stokes equations. In this paper, 
we presented solutions in one-dimensional physical space. 
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1. Introduction 

Gas-liquid flows in small scale geometries have received considerable at- 
tention in the past few years, especially due to the rapid developments in 
micro- and nanofluidics. Recently, experiments have been performed in which 
a segmented flow of gas and liquid is studied in nanochannels [???]. The 
Knudsen number of these flows, i.e. the ratio of the mean free path of the 
particles and the characteristic length, is such that the Boltzmann equation 
needs to be solved to describe the transport processes in the gas phase, while 
usually the incompressible Navier-Stokes equations are sufficient to model 
the liquid phase. While numerical methods for solving the Boltzmann equa- 
tion are well established [? ? ], no efficient methods seem to exist to study 
the flow of a gas at rarefled conditions coupled to the liquid flow described 
by the incompressible Navier-Stokes equations. The aim of this article is to 
take first steps into that direction. 

We consider a liquid droplet inside a gas. A liquid droplet may move, 
for example, due to a pressure difference, so that we need to solve a moving 
interface problem. One may choose different methods to compute these types 
of two-phase flows with moving interface, however, meshfree Lagrangian par- 
ticle methods seem to be one of the preferred choices for such problems. For 
the rarefied gas phase we solve the Boltzmann equation by a DSMC type 
of particle method [? ? ]. In DSMC type of methods the gas particles 
are Lagrangian particles and move with their molecular velocities. How- 
ever, these methods are mesh-based since one divides a computational do- 
main into cells and performs intermolecular collisions of particles inside cells. 
Moreover, macroscopic quantities like fluid density, velocity, temperature, 
etc. are stored in the cell centers. On the other hand, we solve the incom- 
pressible Navier-Stokes equations by the Finite Pointset Method(FPM) [? ? 
], which is a meshfree Lagrangian particle method and is similar in charac- 
ter as Smoothed Particle Hydrodynamics (SPH) [? ]. Within the FPM we 
approximate a liquid domain by Lagrangian particles which move with the 
local fluid velocity. We note that in this article we utilize two types of parti- 
cles which may confuse the readers. Therefore, we call the DSMC particles 
"gas molecules" and the FPM particles "gas- and liquid particles" for the 
gas- and liquid phases, respectively. FPM particles are numerical grid points 
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that move with the fluid velocity and carry all fluid dynamic quantities like 
densities, velocities, pressures, etc. along with them. To solve this type of 
two-phase flow problem one has to decompose the computational domain 
into a liquid and a gas domain. The domain decomposition is performed 
by first generating the entire domain by regular cells which are DSMC cells. 
Then we generate liquid particles overlapping the DSMC cells. The inter- 
face between the two domains is determined from the liquid particles. In 
the one-dimensional case, it is easily determined by identifying the extreme 
ends of the liquid particle positions. For higher- dimensional problems the 
interface is determined by identifying the free surface particles among the 
liquid particles, see [? ? ] for details. 

To the best of the author's knowledge, this is the first attempt to couple 
numerically the incompressible Navier-Stokes equations and the Boltzmann 
equation in the rarefied regime. To validate the coupled solution of the 
equations we have further simulated the gas and liquid phases by solving the 
compressible and the incompressible Navier-Stokes equations, respectively. 
In the latter case the interface conditions are quite standard which have 
been reported in a number of articles in the past few years [????? 
], see also the references therein. We have considered the same examples 
as presented in [? ? ], where the flow has been modeled by the Euler 
or the Navier-Stokes equations. Since we solve both the compressible and 
the incompressible equations by a particle method, we use different flags to 
distinguish the gas and liquid particles. Interfaces between the fluids can 
be determined based on their flags or colors [? ? ]. The flags are defined 
initially, and particles carry them along and leave them unchanged during 
the simulations. 

As will be shown in the article, for small Knudsen numbers the solutions 
obtained from the coupling of the Boltzmann and the incompressible Navier- 
Stokes equations are very close to the solutions obtained from the coupling of 
the compressible and the incompressible Navier-Stokes equations. However, 
the same is no longer true for the larger Knudsen numbers. We present test 
cases with smaller as well as larger Knudsen numbers. 

The paper is organized as follows. In section 2, we introduce the math- 
ematical model for the gas and the liquid phases. In addtion to that we 
derive the interface boundary conditions. In section 3, we describe the par- 
ticle methods for solving the Boltzmann and the Navier-Stokes equations. 
Section 4 is devoted to the coupling algorithms for the gas and the liquid 
phases. The numerical tests are presented in section 5, and some concluding 
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remarks are given in section 6. 



2. Mathematical Model 

2.1. Gas phase: the Boltzmann and the compressible Navier-Stokes equations 

The Boltzmann equation describes the time evolution of a distribution 
function f{t, x, v) for particles of velocity v G at point x G C M^, (s = 
1, 2, 3) and time t G M+. It is given in nondimensional form as 

^ + v-V./ = ^J(/,/), (1) 

where e is the Knudsen number, J(/, /) is the collision operator which is 
given for hard spheres by 

J{fJ)= f I /3(|v-w|,n)[/(v')/(w')-/(v)/(w)]rfnrfw, (2) 

where S"^ is the unit sphere in M^, n G is the unit vector in the impact 
direction, (3 is the collision cross section, /(v ) = fit,:^.,^ ) and analogously 
for /(v) etc. The pairs (v, w) and (v , w ) are the pre- and post colhsional 
velocities of two colliding particles, given by 

V = V — n [n ■ (v — w)] , w = w + n [n ■ (v — w)] . (3) 

The collision operator has five colhsional invariants tp^v) = 1, v, |vp/2 
satisfying 

/ V(v)J(/,/)rfv = 0. (4) 

In other words, J(/, /) locally satisfies the conservation laws for mass, mo- 
mentum and energy. 

The basic quantities of interest are the macroscopic ones, like the density 
p{t, x), mean velocity u = u(t, x) and the total energy E = E{t, x), and are 
defined as 



P 



[ /(^,x,v)rfv, u=i/ v/(t,x,v)rfv (5) 
E=- ^/(t,x,v)dv = -|up + e, (6) 

P JR3 ^ ^ 
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where e is the internal energy, defined by 



e = - f ^^r^/(t,x,v)rfv. (7) 
Moreover, the pressure tensor ip and heat fiux q are defined by 

/ (v-u)®(v-u)/(t,x,v)dv (8) 

q = / ^^^(v-u)/(t,x,v)dv. (9) 

The gas pressure p is defined a.s p = |pe for a monoatomic ideal gas. 
Furthermore, p = pRT holds, where T is the temperature and R is the gas 
constant. For more details we refer to [? ? ]. 

Multiplying the Boltzmann equation by its coUisional invariants and then 
integrating with respect to v over we obtain the following local conserva- 
tion equations 

| + V.(pu) ^ 
^^^"^ + V ■ (pu ® u + (p) = (10) 



dt 
dt 



V ■ [pEu + (/) ■ u + q] = 0. 



For e tending to zero, i. e. for small mean free paths, one can show that 
the phase space distribution function / tends to the local Maxwellian [? ] 

/M(^,x,v) = ^^-^e-^ (11) 

and the system of moment equations fllOl) tends to the compressible Eu- 
ler equations with the closure relations ip = pi and q = 0. This can be 
verified from the asymptotic expansion of / in e, where the zeroth order 
approximation gives the local Maxwellian distribution, and the first order 
approximation [? ] gives the Chapman- Enskog distribution 



/c£;(t,x, v) = /Af(t,x, v) 



q ■ c / |cp 5\ 1 r : c 



hp{RTY\2RT 2) 2 p{RTf 



(12) 
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where c = v — u. At the same time, f llUp tends to the compressible Navier- 
Stokes equations with the closure relations 

if = pi — T, q = — kVT, (13) 

where 



dui duj 2 , 
— - H (V ■ u) 



(14) 



and fi = fi{t, x) and k, = K,{t, x) are the dynamic viscosity and thermal 
conductivity, respectively. They are of the order of e. For example, the first 
approximation for the viscosity and the thermal conductivity of a monatomic 
gas is given by [? ] 



5 mkT 



where k is the Boltzmann constant, and d and m are the mass and the diam- 
eter of the molecules, respectively. In this paper we compute the parameters 
yU, and K from the initial temperature, and use the corresponding values in 
the compressible Navier-Stokes equations. 

Since we solve the Navier-Stokes equations with a meshfree Lagrangian 
particle method, we re-express them in Lagrangian form with respect to the 
primitive variables as 



dpg 

dt 
dt 

dTg 

dt 



+ Va 



[-Pg V 

^vPg 



-V(V 



(16) 



where we have used the index g for the gas quantities, u is the kinematic 
viscosity, d/dt is the material derivative, Cy is the specific heat at constant 
volume, given by |i? for a monoatomic gas. We note that we have expressed 
the internal energy of the gas as eg = CyTg. 

When introducing the specific heat into the energy equation an ideal gas 
was assumed. In addition to the system of equations f|T6|) we consider the 
equation of state 

Pg = PgRTg. (17) 



6 



2.2. Liquid phase: Incompressible Navier- Stokes equations 

We consider incompressible flow inside the liquid phase. The governing 
equations can be obtained from the incompressible Navier-Stokes equations 
by assuming the liquid density pi to be constant. All liquid quantities are 
denoted with the index /. Moreover, one can express the internal energy e; 
for the liquid phase approximately by ci = CpTi, where Cp is the specific heat 
at constant pressure [? ]. We also solve the incompressible Navier-Stokes 
equations by a meshfree Lagrangian particle method. Therefore, we express 
these equations in Lagrangian form 

dt 

dui 
dt 
dT^ 
dt 

where ti is given by ( |T4l) without the divergence of velocity term. In many 
situations, the viscous dissipation term in the energy equation can be ne- 
glected. A detailed discussion about when it is justified to omit that term 
can be found in [? ? ]. In what follows, we limit ourselves to scenarios with 
negligible viscous heating, resulting in a simplified energy equation 

^ = ^V^T,. (19) 

dt Cppi 

2.3. Initial and boundary conditions 

In this paper we consider a one-dimensional computational domain Q = 
[a,b] C M^, where a is always zero and b varies from 10~^ to 10~^. The 
domain is initially decomposed into the gas domain Qg and liquid domain 
Qi = Q\Qg. We consider cases where the liquid domain always remains 
inside of the full domain. Therefore, the boundaries a and b always belong 
to the gas domain. We prescribe boundary conditions for the gas at points a 
and b. Moreover, there are interfaces between the liquid and the gas domains, 
and we have to further specify the interface boundary conditions, described 
in the next subsection. 

In the gas domain Qg we either solve the compressible Navier-Stokes 
equations or the Boltzmann equation. We assume that initially the gas is 







Pi 



(18) 



CpPl 



CpPl 
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in thermal equihbrium with the values Pg{0,x),Ug{0,x) and Tg{0,x), which 
are the initial conditions for the compressible Navier-Stokes equations. If 
we solve the Boltzmann equation in Qg, we prescribe the initial condition 
as a local Maxwelhan with parameters Pg(0, x), Ug(0, x) = {ug{0, x), 0, 0) and 
Tg[0,x). In Qi we solve the incompressible Navier-Stokes equations with 
initial conditions pi{0, x), ui{0, x) and T;(0, x). 

2.4- Interface boundary conditions 

To couple the liquid and gas phases one has to first determine the interface 
between two phases and then prescribe the interface boundary conditions. 
For solving the Boltzmann and the incompressible Navier-Stokes equations, 
we determine the interface as the free surface particles from the liquid do- 
main. In the one dimensional case, we have two interfaces, which are given 
by the liquid particle position at the extreme left (xl) and the liquid parti- 
cle position at the extreme right {xr). They can be tracked at every time 
step. When we solve the transport equations in both phases by the FPM, 
we assign different flags or colors for the particles in the compressible and 
incompressible phases. The particles of each phase carry the color function 
along with them, and the interface can be tracked with the help of the flags 
of the particles. We again determine the interface by identifying the liquid 
particles at the extreme left and right. 

Owing to the kinematic boundary condition at the interface, there is no 
penetration of particles from one phase to the other. This means that the 
convective terms for mass, momentum and energy transport are zero. Hence, 
all fluxes with the multiplicative factors u vanish. Therefore, we have the 
following jump conditions for the momentum and energy fluxes the system 



Here we use a scalar quantity q for the heat flux in the one-dimensional case. 
Moreover, we assume that the velocity and temperature at the interface are 
continuous, i.e. 



m 



{ipuu + q)i 



(<^ll)9 

{ipiiu + q)g. 



(20) 
(21) 



Tl = Tg. 

Then from ([20]), (E]) and ([22]) we have 



Ul = Ug. 



(22) 
(23) 



qi = Qg- 



(24) 
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2.4-1 ■ Interface boundary conditions for the compressible and the incompress- 
ible Navier-Stokes equations 
We use the closure relations f lT3|) and f|T^ in fl20l) and and get the 
following continuity of the fluxes 



Using the divergence-free condition for the liquid we obtain the interface 
boundary condition for the pressure as 



We note that we use the condition ( l23i) together with ( l26i) for the ther- 
modynamic equations. 

2.4-2. Interface boundary conditions between the Boltzmann and the incom- 
pressible Navier-Stokes equations 
An important point to note in relation to the interface conditions be- 
tween the Boltzmann and the Navier-Stokes domain is the fact that while 
the quantities available on the gas side of the interface determine all of the 
quantities needed on the liquid side, the same is not true in the opposite 
direction. Therefore, additional assumptions have to be made when com- 
puting the interfacial phase-space distribution in the gas from the density, 
velocity and temperature in the liquid. This requires a different treatment 
depending on whether information is passed from the gas to the liquid or vice 
versa. As in the previous subsection 12.4.11 we assume the continuity of the 
velocity and the temperature as given by (122|) and (123|) . In this case the clo- 
sure relations ( !T3|) and (fT^ are used only for the fluxes of the incompressible 
Navier-Stokes equations. Therefore, the continuity relations for the fluxes 
differ and are re-expressed in the form 




(25) 



(26) 




(27) 



ij>)l = {Vu)g 



(28) 




(29) 
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where (v^ii)^ and {q)g are computed directly from the moments ([H]) and (jH]), 
respectively, of the solution of the Boltzmann equation. Hence, (12^ . 
and (129|) give the interface conditions from the gas into the liquid, where Tg 
is computed from the solution of the Boltzmann equation. 

On the other hand, the interface boundary condition from the liquid into 
the gas is treated as follows. When gas molecules hit the interface, we apply 
the diffuse reflection condition with thermal accommodation, i.e. particles 
are reflected with the interface temperature and velocity into the gas domain, 
see the detailed descriptions in section 14.21 

3. Numerical methods 

3.1. Particle method for the Boltzmann equation 

For solving the Boltzmann equation we have used a variant of the DSMC 
method [? ], developed in [? ? ]. The method is based on the time split- 
ting of the Boltzmann equation. Introducing fractional steps one first solves 
the free transport equation (the coUisionless Boltzmann equation) for one 
time step. During the free flow, boundary and interface conditions are taken 
into account. In a second step (the collision step) the spatially homoge- 
nous Boltzmann equation without the transport term is solved. To solve the 
homogeneous Boltzmann equation the key point is to find an efficient par- 
ticle approximation of the product distribution functions in the Boltzmann 
collision operator given only an approximation of the distribution function 
itself. To simulate this equation by a particle method an explicit Euler step 
is performed. To guarantee positivity of the distribution function during the 
collision step a restriction of the time step proportional to the Knudsen num- 
ber is needed. That means that the method becomes exceedingly expensive 
for small Knudsen numbers. 

3.2. FPM for the compressible Navier-Stokes equations 

We solve the Navier-Stokes equations (fT6|) by the FPM. As already pointed 
out, the FPM is a meshfree Lagrangian particle method, where we approx- 
imate the spatial derivatives with the help of the weighted least squares 
method. In order to solve these equations by FPM, one first fills the compu- 
tational domain by particles which can be considered as moving numerical 
grid points and then approximates the spatial derivatives occurring on the 
right hand side of f|T6|) at each particle position from its neighboring par- 
ticles. This reduces the system of partial differential equations f lTB]) to a 
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system of ordinary differential equations with respect to time. We solve the 
resulting ODE system with the help of a two-step Runge-Kutta method. 
The time steps for the compressible as well as the incompressible Navier- 
Stokes equations are restricted by the CFL condition and by the value of the 
transport coefficient max [fig/pg, Hg/{pgCy), ni/{piCp)]. We refer to our earlier 
papers [? ? ] for the details of the least squares approximation of the spa- 
tial derivatives. We note that we need to introduce a particle management 
scheme during the simulations. Because of the Lagrangian description of 
the method particles may accumulate or may thin out causing holes in the 
computational domain. This gives rise to some instabilities of the method. 
Therefore, we have to add or remove particles. In the one-dimensional case 
this task is quite simple. If the distance between a particle and its nearest 
neighbor is larger than 1.2 times the initial spacing of particles, we add a 
new particle in the center. On the other hand, if two particles are closer 
than 0.2 times the initial spacing we remove both of them and add a new 
particle at the mid point. However, these adding and remove factors may 
depend upon users. The fluid dynamic quantities of newly added particles 
are approximated from their neighbor particles with the help of the least 
squares method. 

3.3. FPM for the incompressible Navier-Stokes equations 

The incompressible Navier-Stokes equations (fT5]l are also solved by the 
FPM. Since we consider one-dimensional physical space, the solution of the 
incompressible Navier-Stokes equations is simple. For higher-dimensional 
physical spaces we refer to [? ? ]. 

Taking the partial derivatives with respect to x on both sides of the 
momentum equation flTS]) and using the divergence- free constraint, we obtain 
the one-dimensional Laplace equation for the pressure 

with the Dirichlet boundary conditions pi and on the left and right in- 
terfaces xl and xr, respectively. 

Suppose we have N liquid particles x/i, ■ ■ ■ ,xin in fli. We have to solve 
the corresponding transport equations for every liquid particle. Let xu be 
the position of an arbitrary particle. The new pressure at time level t""*"^ = 
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(n + l)At, where At is the time step, can be computed exphcitely by 



Pli ~ ^n-^li^ ' y"^^) 



^R -^L ^R 



The new velocity is given by 

Pl -^R -^L 

The energy equation for the hquid is solved in the same way as in the 
case of the gas. 

We note that the equation for the velocity does not contain any contri- 
butions from the viscous term, which is due to the one- dimensional nature of 
the problem. The thermal diffusivity Ki/{piCp) of the liquid is much smaller 
than that of the gas. Therefore, the time step for the gas phase is smaller 
than the one for the liquid. So, we use an explicit Euler step for the time in- 
tegration to solve the energy equation of the liquid phase. However, the time 
steps are chosen as the smallest of the values obtained for the two phases. 

Finally, we compute the new liquid particle positions by 

= xi^ + Atul. (33) 



4. Coupling of the gas and the liquid phase 

The model for the coupling of the two phases to be presented in the fol- 
lowing contains some essential features of gas-surface interactions, but also 
relies on some simplifications. Essentially, we neglect evaporation and con- 
densation phenomena, i.e. we assume mass conservation for both the gas 
and the liquid. While this could be an unsupported assumption for cases 
in which rarefaction is due to the reduced density of the gas phase, in the 
cases we study the reduction of the length scale down to micrometer dimen- 
sions causes rarefaction. Still, the results to be presentented in the following 
sections show that there exist local variations of the gas density and temper- 
ature. Therefore, strictly speaking, while the global (integrated) evaporation 
and condenstion mass fluxes could be very small, they may be of relevance 
locally. Nevertheless, two arguments show that even our comparatively sim- 
ple model for gas-liquid interaction can make practically relevant predictions 
about gas-liquid flows. First of all, we study fast processes occuring on time 
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scales below on microsecond where phase-change phenomena are usually neg- 
hgible. Moreover, there exist non-volatile liquids for which even on longer 
time scales evaporation does not play a role [? ]. Therefore, with the neces- 
sary precations taken, the model presented in the following should provide a 
realistic picture of a specific class of gas-liquid flows. 

4.1. Coupling of the compressible and incompressible Navier-Stokes equations 

We initially generate liquid and gas particles in the computational do- 
main. We assign all initial values such as density, velocity, temperature, etc. 
to all particles. Moreover, we also assign different flags or colors to liquid and 
gas particles. For instance, we define a flag equal to 1 for liquid particles and 
2 for gas particles. These flags remain unchanged throughout the simulation. 
The gas equations are solved by gas particles and the liquid equations by liq- 
uid particles. We approximate the spatial derivatives at a gas particle from 
its neighbors and exclude all liquid particles from the neighbor list. This 
situation occurs near the interface. Similarly, we exclude the gas particles 
from the neighbor list of liquid particles. The interface particles at xl and 
Xfi are determined after every time steps via the leftmost and rightmost po- 
sitions among the liquid particles. For higher dimensional cases, we need to 
spend more effort to determine the interface particles, but this is also quite 
straightforward, see [? ? ]. 

To couple both types of equations, we first solve the compressible Navier- 
Stokes equations for the gas particles. We consider the velocity ui and the 
temperature Ti of the interface particles as Dirichlet boundary conditions. 
This gives the coupling conditions from the liquid phase into the gas phase. 

The coupling from the gas phase into the liquid phase is obtained as fol- 
lows. After computing new quantities for the gas particles, we approximate 
the pressure and temperature at the interface positions xl and x^. To ap- 
proximate the pressure according to (!27|) . we first extrapolate the pressure 
from the neighboring gas particles using the least squares method. Then 
we subtract the term containing the derivative of Ug which is also computed 
from the neighboring gas particles from the least squares method. Subse- 
quently we obtain the liquid pressure from f l3ip and the velocity from f l32p 
after determining the pressure on the interface. 

The approximation of the temperature at the interface is slightly different. 
At the interface the temperature and the heat flux are continuous according 
to f l23|) and fl26|) . implying a jump condition for the temperature gradient. We 
apply a two-sided interpolation satisfying these jump conditions as follows. 
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Let xj be a interface particle with the temperature T/. The particle x/ has m 



neighboring gas particles Xgi,Xg2, 



,x 



gm 



and n neighboring liquid particles 



Xii,xi2, ■ ■ ■ , xin- We write the Taylor expansions of Tu and T„j around T/ as 



T • 

-^gi 



Tu 



dT 

T. + dxA^) 

dTi- 
dx, 



+ e 



gii 



m 



(34) 
(35) 



where Cgi and en are the errors in the Taylor's expansion at the interface point 
xj, dxii = xii — xj, and dxgi = Xgi — xj. The continuity of the temperature 
field {Ti = Tg = Tj) is reflected in the first terms of the right-hand side 
of flMl) and fl35l) . In addition to that, we add the jump condition fl26l) in 
the system of equations (134|) and ( 135|) as a constraint. Then we obtain the 
interface temperature T/ after solving the m + m + 1 system of equations from 
minimizing the errors [? ? ]. The system of equations can be re-expressed 
in the following matrix form 



/ 1 dxgi \ 



1 dxg^ 



1 



1 

\0 





Hi 



dxii 

dXin 



dTg_ 

dXa 



\ ^ / 

\ dxi I 



J 



( Tgl \ 



T 

gm 

Tn 



'''In 

\ / 



+ 



'-'gm 

en 



ein 
\ ei J 



(36) 



Here the cj denotes the error in the continuity of of the heat flux at the inter- 
face. This constraint interpolation gives accurate solutions of the diffusion 
equation with discontinuous diffusion coefficients [? ]. After determining 
the temperature at the liquid interface, we use them as Dirichlet boundary 
conditions for the thermodynamic equation (fT9!) and obtain the temperature 
for the liquid phase. 

This provides the coupling conditions from the gas phase into the liquid 
phase. To summarize the above descriptions, we present the following al- 
gorithm for coupling the compressible and the incompressible Navier-Stokes 
equations. 

4.1.1. Coupling Algorithm I 

(i) Generate initial particles with flags as liquid and gas and prescribe the 
initial values. 
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(ii) Determine the interface particles xl and xr from the set of hquid parti- 
cles. 

(iii) Solve the compressible Navier-Stokes equations with the Dirichlet bound- 
ary conditions at the interface (liquid) particles. 

(iv) Interpolate the pressure and the temperature to the interface particles. 

(v) Solve the incompressible Navier-Stokes equations. 

(vi) Move all particles with their velocities. 

(vii) Add or remove particles if necessary, see subsection 13.21 

(viii) Goto (ii) and repeat until the final time is reached. 

4-2. Coupling of the Boltzmann and the incompressible Navier-Stokes equa- 
tions 

The coupling procedure for the Boltzmann and the incompressible Navier- 
Stokes equations is different from the previous one. The Boltzmann equa- 
tion is a mesh-based method since gas molecules have to be sorted in cells 
for the intermolecular collisions. Moreover, we sample and store the macro- 
scopic quantities at the cell centers. As already described, the incompressible 
Navier-Stokes equations are solved by a meshfree method. So, we need to 
couple the mesh-based and the meshfree method. In this coupling proce- 
dure, we first generate regular cells in the entire domain. In addition to that, 
we fill the liquid domain by liquid particles. The liquid particles overlap 
with the Boltzmann cells. We generate gas molecules outside of the liquid 
domain. Then we initialize the velocity distribution as a Maxwellian with 
the initial parameters pg{0,x),Ug{0,x),Tg{0,x) in all Boltzmann cells. As a 
consequence we have cells with and without gas molecules. If a cell contains 
no gas molecule we deactivate it, while cells with gas molecules are active 
cells. We consider a fixed number of cells for the Boltzmann solver. If the 
cell size is larger than the mean free path, we perform a refinement such that 
the new cell size is smaller than the mean free path, see [? ] for details. 

Moreover, we prescribe initial conditions for the liquid particles. As in 
the earlier case we can determine the left and right boundaries of the liquid 
domain. First, we solve the Boltzmann equation. After the free flow step 
in the Boltzmann solver we apply the boundary conditions. If gas molecules 
cross the interface to the liquid, they lose the memory of their velocity and 
are reflected diffusively with a Maxwellian distribution having the local ve- 
locity and temperature parameters of the interface particles. This gives the 
coupling conditions from liquid phase into the gas phase. 
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Similarly, if the domain boundaries at a and b are considered as solid 
walls and gas molecules cross these boundaries, they also lose the memory 
of their velocity and reflect diffusively with a Maxwellian distribution having 
the wall temperature and velocity parameters. 

The coupling procedure from the gas phase into the liquid phase is similar 
to the coupling of the compressible and incompressible Navier-Stokes equa- 
tions. At the end of a Boltzmann solver step we compute the stress tensor 
and heat flux from the gas molecules and store them at the cell centers. 
Moreover, we also compute the density, velocity, temperature and pressure 
and store them at the cell centers. Then we approximate the pressure at the 
interface from the neighboring activated cell centers of the Boltzmann solver 
using equation f l28p . The temperature at the interface is approximated sim- 

with the conditions (ES]) and f l2^ . 



ilarly as in subsection 14. ![ with the conditions f l23|) and f l29|) . Note that the 
system of equations differs slightly from f p6|) in this case and is re-expressed 
as 
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It is well known that in all DSMC type solvers there are some statistical 
fluctuations in the solutions. If we pass these fluctuating data, they desta- 
bilize the Navier-Stokes solver. Therefore, we need a smoothing operator. 
Here we have used the Shepard interpolation. For example, for a function i/) 
at a cell center x, the Shepard interpolation is defined as 
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where m is the number of neighboring active cell centers Xj, and Wi are the 
weight functions depending on the distance between x and Xj. We have 
chosen a Gaussian weight function given by 
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with h being the radius of interaction which is three times the cell distance, 
a is a positive constant, taken as 2, which, however, can be varied. We 
have already used this type of smoothing for coupling the Boltzmann and 
the Navier-Stokes equations [? ], with the result that for small Knudsen 
numbers the locally smoothed Boltzmann solutions are compatible with the 
non-smoothed Navier-Stokes solutions. This means that there is negligible 
smearing of the solutions with this type of local smoothing. 

After approximating the pressure and temperature at the gas-liquid inter- 
face, we solve the incompressible Navier-Stokes equations in a similar manner 
as in subsection 14.11 This completes the description of the coupling from the 
Boltzmann domain into the incompressible Navier-Stokes domain. 

Recall that in the liquid phase a Lagrangian particle method is applied. 
Therefore, after moving the liquid particles, we may have some gas molecules 
inside the liquid domain. Since we switch to the Boltzmann solver after a 
Navier-Stokes step has been completed, we do not reflect them from liq- 
uid interface at this level, but we treat them after the free-flow step of the 
Boltzmann solver. 

Summarizing the above we present the following algorithm for the cou- 
pling of the Boltzmann and the incompressible Navier-Stokes equations. 

4.2.1. Coupling Algorithm II 

(i) Generate Boltzmann cells in the entire domain and generate liquid parti- 
cles which overlap the Boltzmann cells. 

(ii) Generate gas molecules outside the liquid domain according to a Maxwellian 
distribution with the initial parameters and prescribe the initial conditions 
for the liquid particles. 

(iii) Determine the interface particles xl and xr for the liquid phase. 

(iv) Solve the Boltzmann equation with the gas-liquid interface taking the 
role of a moving wall. 

(v) Compute the moments in the Boltzmann cells and smooth them accord- 
ing to dil). 

(vi) Interpolate the pressure and temperature to the interface particles. 

(vii) Solve the incompressible Navier-Stokes equations. 

(viii) Move the liquid particles with their velocities. 

(ix) Add or remove liquid particles, if necessary, see subsection 13.21 

(x) Goto (iii) and repeat until the final time is reached. 
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5. Numerical Tests 

As we have already mentioned, we consider the computational domain 
Q = [a,b], where a = in all cases and b varies. In Qi we solve the 
incompressible Navier-Stokes equations whereas in Qg we solve either the 
Boltzmann equation or the compressible Navier-Stokes equations. Then we 
compare the results obtained from both coupling algorithms. We consider 
three different test cases. In all cases a hquid droplet is surrounded by a com- 
pressible gas inside the computational domain, where the boundary points a 
and b always belong to Qg. While coupling the compressible and incompress- 
ible Navier-Stokes equations we consider a total number of initial particles 
N = 200 for both liquid and gas domains. The corresponding initial grid 
spacing is Ax = 1/N. To determine the neighboring particles in the FPM, 
the interaction radius h is taken equal to 3 x Ax. 

For the coupling of the Boltzmann and the incompressible Navier-Stokes 
equations we generate = 200 regular cells in Q and add liquid particles 
which cover the liquid domain Qi. Note that in algorithm I, Qg and Qi are 
disjoint, but same is not true in algorithm II. The initial grid spacing of 
liquid particles is again equal to The Boltzmann cell size is refined 

according to the size of the mean free path, see [? ] for details. As we have 
mentioned initially, the hard-sphere model for the collision cross section is 
employed for the Boltzmann equation. Moreover, we have used the following 
parameters for all test cases. The gas is Argon with a molecular mass m = 
6.63 X 10"^^ kg. For the Boltzmann constant we have k = 1.38 x 10~^^ J/K, 
for the molecular diameter d = 3.68 x 10~^° m, and we obtain the gas constant 
R = 208 J/{kgK). The dynamic viscosity and thermal conductivity for the 
compressible Navier-Stokes equations are assumed to be constant and are 
evaluated with the initial temperature according to ( 1T5|1 . 

For the liquid phase, in all three cases we assume an initial temperature 
equal to 298 K. The thermal conductivity and the specific heat capacity of 
the liquid are taken to be those of water, giving values of ki = 0.63 J/ {mKs) 
and Cp = 4181 J/{kgK), respectively. 

5.1. Test 1 

We first study a liquid droplet with nonzero initial velocity. This test 
case has been considered by Caiden et al [? ] for the inviscid case. We first 
consider the interval Q = [0 m, 1 x 10~^ m]. Initially, a liquid droplet occupies 
the domain Qi = [4 x 10~^ m, 6 x 10"^ m], while the gas occupies the rest 
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of the domain. The gas has the initial values Pg = 1.226 kg/m^,Ug = m/s 
and Pg = 1 X 10^ Pa. The temperature is obtained from the equation of 
state. The initial values for the liquid phase are pi = 1000 kg/m^,ui = 
100 m/s,Ti = 298 K. The initial pressure in the liquid is taken to be the 
same as the initial pressure in the gas, however, this will change due to the 
coupling of the phases. 

The boundary values for the compressible Navier-Stokes equations are 
given by Ug(t,a) = Ug(t,b) = Ug{0,x) and Tg(t,a) = Tg(t,b) = Tg{0,x). 
Moreover, for the Boltzmann solver, if gas molecules cross the physical do- 
main, we reflect them back in the same was as in the case of the gas-liquid 
interface boundary conditions, using the boundary velocity and temperature 
Ug{t, a) = Ug(t, b) = Ug{0, x) and Tg(t, a) = Tg{t, b) = Tg(0, x). 

The characteristic length is taken to be the length of the droplet which 
is equal to 2 x 10~^ m. The corresponding Knudsen number is e = 0.0045. 
The simulations are stopped after t = 5.2 x 10~^ s. 

The initial number of gas molecules is 5000 per cell of size Ax for the sim- 
ulation of the Boltzmann equation. We observe in figure [T] that the solutions 
obtained from the coupling of the Boltzmann and the incompressible Navier- 
Stokes equations (Algorithm II) match with the solutions obtained from the 
coupling of the compressible and incompressible Navier-Stokes equations (Al- 
gorithm I). 

As in [? ] we notice that the liquid droplet moves to the right causing 
a compression wave in the gas ahead of it and an expansion wave behind it. 
Moreover, the gas is heated ahead of droplet and cooled behind it. To a very 
good approximation, the pressure profile inside the droplet is linear. This 
is expected and can be explained by the following argument. The liquid is 
decelerated by compressing the gas volume at its right. By the equivalence 
principle of general relativity, an acceleration (or deceleration) is equivalent 
to the action of a gravitational field. Therefore a linear hydrostatic pressure 
profile as in a quiescent liquid in a gravitational field is found. 

Next, we consider a 10 times smaller domain = [0, 1 x 10"^ m] and the 
same initial values as above. Initially a liquid droplet occupies the domain 
fii = [Ax 10~^ m, 6 X 10~^ m], while the gas occupies the rest of the domain. 
The characteristic length is taken as the size of the droplet, which is equal to 
2 X 10~^, giving a Knudsen number e = 0.045. The simulation is stopped after 
a time t = 5.2 x 10~^ s. The results are plotted in figure |2] and show similar 
trends as the curves in figure [1] However, at such values of the Knudsen 
number first deviations from the results of continuum theory become visible. 
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Figure 1: Testl: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a Knudsen number of 0.0045. The solid lines are results 
from Algorithm I, the circles are from Algorithm II. 
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Figure 2: Testl: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a Knudsen number of 0.045. The solid lines are results 
from Algorithm I, the circles are from Algorithm II. 
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Figure 3: Testl: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a Knudsen number of 0.45. The solid lines are results 
from Algorithm I, the circles are from Algorithm II. 
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We have further increased the Knudsen number by decreasing the com- 
putational domain and characteristic length by a factor of 10. In this case 
we have e = 0.45, and simulation was stopped after t = 5.2 x 10"^° s. We 
have plotted the results in figure [31 For this Knudsen number we see that the 
solutions from the coupling algorithm I significantly deviate from the ones 
obtained from coupling algorithm II. This is, as expected, due to the failure 
of the compressible Navier-Stokes equations for larger Knudsen numbers. A 
constant velocity and a linear pressure profile are obtained inside the liq- 
uid. It also becomes apparent that the solution of the Boltzmann equation 
yields much more pronounced jumps in the temperature fields than that of 
the Navier-Stokes equations, similar as in the case of Sod's ID shock tube 
problem [? ]. 

5.2. Test 2 

The second test case is also taken from paper of Caiden et al [? ] , where 
the authors have considered only the Euler equations without the energy 
equation. We again consider the computational domain Q = [0 m, 1 x 
10~^ m]. The liquid initially occupies the domain Qi = [A x 10~^ m, 6 x 
10"^ m], and the rest of the domain is filled with gas. A shock wave is 
initially located at x = 1 x 10~^ m, with a post shock state pg{0,x) = 
148407.3 Pa, Ug{0, x) = 89.981 m/s and pg{0, x) = 2.214 kg/m^. On the right 
of X = 1 X 10"^ m the initial state is given as Pg{0, x) = 98066.5 Pa, Ug{0, x) = 
m/s and pg{0,x) = 1.58317 kg/rn?. The initial temperature of the gas 
is again computed from the equation of state. The characteristic length 
is equal to 2 x 10~^ m. The pre- and post-shock Knudsen numbers are 
0.00259 and 0.00696, respectively. The initial conditions for the liquid are 
Pi = 98066.5 Pa,ui = m/s,Ti = 298 K. At first, we first consider a liquid 
density of pi = 1000 kg/m?. 

The boundary conditions for the compressible Navier-Stokes equations 
are Ug(t,a) = 89.981 m/s,Tg(t,a) = T(0,a), Tg{t,h) = Tg{0,b), and we 
assume a vanishing velocity gradient at x = b. 

For the Boltzmann equation we generate gas molecules in two ghost cells 
at the left boundary x = a = according to a Maxwellian with the initial 
parameters pi{0,a),Ug{0,a) = (89.981,0,0) m/s and Tg{0,a). Similarly, we 
generate gas molecules in two cells at the right boundary x = b according to 
a Maxwellian distribution with the parameters pg{0,b),Tg{0,b). The mean 
velocity Ug{t,b) is extrapolated from the neighboring cells at the left. After 
the free flow step, if gas molecules leave the boundaries at a and b we delete 
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Figure 4: Test 2: Logarithm of density (top left), pressure(top right), velocity (bottom 
left) and temperature (bottom right) for a liquid density pi ~ 1000 kg/vT" at time t = 
1.75 X 10^^ s. The solid lines are results from Algorithm I, the circles are from Algorithm 
IL 
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them. We note that we stopped the simulation before the wave reached the 
right boundary so that the mean velocity is still zero at this time. We have 
plotted the results at time t = 1.75 x 10~^ s in figure HI One can observe that 
the shock wave travels to the right hitting the droplet, causing refiected as 
well as transmitted waves. The refiected wave is clearly visible in the pressure 
and the temperature plots of figure HI while both of the waves can be seen 
in figure O Again, the pressure profile is linear inside the droplet, however, 
with a reversed slope compared to the first test case. This originates from 
the fact that due to the interaction with the shock wave the accelaration of 
the droplet is now positive. The solutions obtained from both algorithms are 
very close to each other. 

In figure we have plotted the results from the same simulations with 
an liquid density pi = 10 kg/m^ at the same time t = 1.75 x 10"'' s. It 
becomes apparent that for this lower liquid density the droplet has a lower 
inertia and is more easily displaced to the right. Compared to the initial 
state, the physical fields of the fiuids have now been altered in almost the 
entire domain. Again, we see a good agreement of the solutions obtained 
from the coupling algorithms I and II. 

The above results are for the continuum gas regime. As expected, for the 
rarefied regime we see a discrepancy of the solutions obtained from algorithms 
I and II as in Test 1 for the largest Knudsen number. We do not explicitly 
present the corresponding results in this article. 

5.3. Test 3 

In the final test case we again consider the same computational do- 
main as in Test 2. The liquid droplet is initially at rest, occupying the 
domain Qi = [2 x 10~^ m, 3 x 10~^ m], while the rest of the domain is 
filled with gas. A similar test case has been studied in [? ] for a unit 
interval. Initially, the gas is in thermal equilibrium with the initial state 
Ug{0,x) = (0,0,0) m/s, Tg{0,x) = 298 K, and a spatially varying density 
Pg{0,x). We assume pg{0,x) = 1 kg/rn? at the left of the droplet and a four 
times smaller value at the right of it. This means that the pressure at the 
left of the droplet is four times larger than at the right. The characteristic 
length is 1 X 10~^ m, and the Knudsen numbers at the left and right of the 
droplet are 0.011 and 0.044, respectively. The initial values for the liquid are 
u/(0,x) = Om/s, Ti(0,x) = 298 K, and the pressure is equal to the initial 
value in the gas at the right side of the droplet. We set the liquid density 
pi = 10 kg/m?. The boundary conditions are the same as in Test 1. In 
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Figure 5: Test 2: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a liquid density pi ~ 10 kg/m^ at time t = 1.75 x 10~^ s. 
The solid lines are results from Algorithm I, the circles are from Algorithm II. 
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Figure 6: Test 3: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a liquid density pi = 10 kg/m^ at time t — 2.218 X 
10~^ s. The solid lines are results from Algorithm I, the circles are from Algorithm II. 
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figure El we have plotted the density, pressure, velocity and temperature after 
a short time t = 2.218 x 10~^ s. We see that the droplet has been slightly 
pushed to the right. The pressure on the left is higher. The gas has started 
to cool down on the left, and on the right it is getting compressed, resulting 
in a temperature increase. 

At a later time t = 5.678 x 10~^ s (figure [7]) we observe that the compres- 
sion on the right has further increased and the gas has been further heated. 
Now the pressure on the right is higher causing a slowing down of the droplet. 
Going along with the transition from droplet accelaration to deceleration is 
a reversal of the slope of the pressure profile inside the droplet. Later it will 
bounce back to the left. We observe that at both time levels the solutions 
obtained with both algorithms are in good agreement. 

At a still later time [t = 9.938 x 10"^ s) we see that the droplet is already 
pushed to the left, as displayed in figure El At this time we see some discrep- 
ancy between algorithms I and II. This may be due to statistical fluctuation 
in the Boltzmann domain. However, the solutions are still comparable. If 
we compute further, the droplet will bounce back at the right boundary, and 
this oscillatory process will continue. 

Similarly, for the rarefied gas regime, the solutions from Algorithms I and 
II do not match due to the breakdown of the continuum hypothesis inside the 
gas domain. As in Test 2 we do not present the results for larger Knudsen 
numbers. 

6. Conclusion and Outlook 

We have presented a coupling algorithm for a moving liquid phase inside 
a rarefied gas. The liquid phase is modeled by the incompressible Navier- 
Stokes equations, while the rarefied gas phase is modeled by the Boltzmann 
equation. The transport equations in the liquid are solved by a meshfree 
Lagrangian particle method, the Boltzmann equation by a DSMC type of 
particle method. Liquid particles overlap with the Boltzmann cells and the 
gas-liquid interface is determined from the particles at the boundary of the 
liquid domain. Interface boundary conditions are derived and a coupling 
algorithm for solving the Boltzmann equation in combination with the in- 
compressible Navier-Stokes equations is presented. To validate the coupling 
between the rarefied gas and the liquid we have modeled the gas by the com- 
pressible Navier-Stokes equations. Both the compressible and incompress- 
ible Navier-Stokes equations are solved by a meshfree Lagrangian particle 
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Figure 7: Test 3: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a liquid density pi = 10 kgjwi? at time t = 5.678 X 
10~^ s. The solid lines are results from Algorithm I, the circles are from Algorithm II. 
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Figure 8: Test 3: Logarithm of density (top left), pressure(top right), velocity (bottom left) 
and temperature (bottom right) for a liquid density pi = 10 kgjwi? at time t — 9.938 X 
10~^ s. The solid lines are results from Algorithm I, the circles are from Algorithm II. 
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method, called the Finite Pointset Method (FPM), where the gas and liquid 
particles are distinguished by assigning different ffags. We have also derived 
and implemented the interface boundary conditions for the compressible and 
incompressible Navier-Stokes equations which are well understood and widely 
used. In the continuum regime we have shown that the coupled solutions of 
the compressible and incompressible Navier-Stokes equations match well with 
the coupled solutions of the Boltzmann and the incompressible Navier-Stokes 
equations. 

Future work will concentrate on the extension of the code to higher- 
dimensional problems, especially to the simulation of nano bubbles sur- 
rounded by an aqueos phase. Moreover, the method will be extended to 
simulate two-phase flows with evaporation and condensation [? ? ]. 
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